#include "Multigrid_Solver.h"
#include "Multigrid_Solver_1D.h"
#include "Multigrid_Solver_2D.h"
#include "Multigrid_Factory.h"

double D1basic_f2rd(double * P)
{
    return  -(sin(P[0])-cos(P[0])*cos(P[0]))*exp(sin(P[0]));
};

double D1basic_f(double *P)
{
    return exp(sin(P[0]));
};

double D2basic_f2rd(double * P)
{
    return  (-2*sin(P[0])*sin(P[1])+sin(P[1])*cos(P[0])*cos(P[0])+sin(P[0])*cos(P[1])*cos(P[1]))*exp(sin(P[0])*sin(P[1]));
};

double D2basic_f(double *P)
{
    return exp(sin(P[0]*sin(P[1])));
};

int main()
{
    Factory1DMGSolver Solver1D;
    Factory2DMGSolver Solver2D;
    
    FactoryMGSolverBase *_p_fac_base = &Solver1D;
    Multigrid_Solver * _arena1D = _p_fac_base ->CreateMGSolver();
    _p_fac_base = &Solver2D;
    Multigrid_Solver * _arena2D = _p_fac_base->CreateMGSolver();

double _a[1]={0};
int n=128;
    VectorXd u_,f_,r_,r_origin,acu;
    acu.setZero(n+1);
for(int i=0;i<=n;i++)
{   
    _a[0]=1.0*i/n;
    acu(i)=D1basic_f(_a);
};

int _Rtype;
int _Itype;
std::cout<<"Give me n and Rtype and Itype,:"<<std::endl;
std::cin>>n;
std::cin>>_Rtype;
std::cin>>_Itype;
_arena1D->Initialize(n,D1basic_f2rd,D1basic_f);
_arena1D->setRestircType(_Rtype);
_arena1D->setIntpType(_Itype);
_arena1D->setSmoothType(2);

double acu_norm=acu.lpNorm<Eigen::Infinity>();
double h;

u_=_arena1D->get_u_origin();
f_=_arena1D->get_f_origin();
h = _arena1D->get_h();
r_origin=_arena1D->ResidualSolver(u_,f_,h);
double r_origin_norm=r_origin.lpNorm<Eigen::Infinity>();
int cnt;
for (cnt=1;cnt<=100;cnt++)
{
    u_=_arena1D->Vcycle(u_,f_,h);
    r_=_arena1D->ResidualSolver(u_,f_,h);   
    if (r_.lpNorm<Eigen::Infinity>() <= r_origin_norm*1e-8)
    break;
};

std::cout<<"n= "<<n<<std::endl;
 std::cout<<"cnt= "<<cnt<<std::endl;
//  std::cout<<wucha<<std::endl;
std::cout<<"Residual norm_inf="<<r_.lpNorm<Eigen::Infinity>()<<std::endl;

std::cout<<"relative accu= "<<r_.lpNorm<Eigen::Infinity>() /r_origin_norm<<std::endl;
    

};

